#PSD on profilometer of slumped parts
%reset
%load_ext autoreload
%autoreload 2
%matplotlib qt
#plot psd of calibration flat to get signature of instrument
#compare profilometer data of OAB103 different measurement with different coatings
#add CCI OAB mandrel data from first set
#WFS data on profilometer and mid-frequencies
#compare psd of slumping models at different temperatures
%qtconsole
import matplotlib.pyplot as plt
import numpy as np
import os
from astropy.io import fits
from dataIO.span import span
from pyGeneralRoutines.fn_add_subfix import fn_add_subfix
from IPython.display import display
#from skimage import data as skidata
#from skimage import img_as_float
#from plot_grains import *
#from radialProfiles import *
#from peak_finder import *
from pyProfile.psd import plot_sig_psd
from psd2d import *
from pySurf.points import *
## trick to import from a non released library
import sys
sys.path.append(r'..\OP1surface')
from fit_cylinder import *
outfolder=r'PSDtest'
#this is the simplest way to make ipython print numbers in non scientific format
%precision '%f'
np.set_printoptions(formatter={'float_kind':lambda x: '%f'%x})
def correctedPSD(pts,crop=None,step=None,nsigma=None,outname=None):
"""Perform a set of operations on pts and plot intermediate
steps:
rebin based on step, this must be provided to convert from pointcloud to matrix
remove outliers line by line on the base of sigma
if outname is provided, save figures.
Return 1d frequency and matrices of psds along axis before and after correction.
"""
if outname is not None:
outfolder=os.path.dirname(outname)
outname=os.path.basename(outname)
plt.figure(1)
plt.clf()
plot_points(pts,scatter=1,aspect='auto')
display(plt.gcf())
if crop is not None:
pts=crop_points(pts,xrange=crop[0],yrange=crop[1])
pts=level_points(pts)
print "rebin with step %s,%s"%step
# rebin
## create a grid for bins
s=span(pts[:,0:2],axis=0)
xb=np.arange(s[0][0],s[0][1],step[0]) #x and y bins
yb=np.arange(s[1][0],s[1][1],step[1])
p,x,y=rebin_points(pts,bins=(xb,yb),matrix=True)
plt.figure(2)
plt.clf()
plt.imshow(p,extent=np.hstack([span(x),span(y)]),aspect='equal',
vmin=-10,interpolation='none',
origin='lower')
plt.colorbar()
if outname is not None:
plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,
'_data','.png')))
display(plt.gcf())
if nsigma is not None:
#remove outliers and plot comparison
pl=remove_outliers2d(p,fignum=3,nsigma=nsigma)
if outname is not None:
plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,'_datacorr','.png')))
display(plt.gcf())
else:
pl=p
f,pp1,pp2=compare_2dpsd(p,pl,x,y,fignum=5)
if outname is not None:
plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,'_2dPSDcorr','.png')))
display(plt.gcf())
ap1=avgpsd2d(pp1)
ap2=avgpsd2d(pp2)
plt.figure(4)
plt.clf()
#plt.title('Calibration mandrel')
plt.plot(f[1:],ap1[1:-1],'r',label='Raw')
plt.plot(f[1:],ap2[1:-1],'b',label='Outlier removed')
plt.loglog()
plt.ylabel('um^2')
plt.xlabel('mm^-1')
plt.legend(loc=0)
plt.grid(1)
if outname is not None:
plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,'_PSD','.png')))
display(plt.gcf())
return f,ap1,ap2
#try to repeat results with a function
df=r'../OP1surface/measure_data/calibration_m3/2013_12_23/OAB150x150_01_Height.txt'
pts=get_points(df,center=(0,0),scale=(-1,1.,1))
crop=[[-60,60],[-60,60]]
step=(0.5,0.5)
nsigma=1.5
#pts=get_points(df,center=(0,0),scale=(1,1,1),delimiter=' ')
f,ap1,ap2=correctedPSD(pts,crop=crop,step=step,nsigma=nsigma)
plt.savefig(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.png')))
np.savetxt(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.dat')),
np.vstack([f[:],ap1[:-1],ap2[:-1]]).T)
#try to repeat results with a function
df=r'C:\Users\Vincenzo\Google Drive\POOL\OP1surface\OAB1\xysurface150x150_Height_transformed_deltaR_4inch.dat'
crop=None
step=(0.5,0.5)
nsigma=1.5
pts=get_points(df,center=(0,0),scale=(1,1,1),delimiter=' ')
f,ap1,ap2=correctedPSD(pts,crop=crop,step=step,nsigma=nsigma)
plt.savefig(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.png')))
np.savetxt(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.dat')),
np.vstack([f[:],ap1[:-1],ap2[:-1]]).T)
#Here I start from raw data, need to fit cylinder
df=r'..\OP1surface\OP1\09_op1_yxsurf_Height_transformed.dat'
#measure_data\2015_07_10\06_OP1_xyscan_Height.txt'
crop=[[-23,23],[-50,50]]
step=(0.5,0.5)
nsigma=1.5
pts=get_points(df,center=(0,0),scale=(1,1,1),delimiter=' ') #z-1 per mandreino
pts=crop_points(pts,crop[0],crop[1])
plt.clf()
a=plot_points(pts,scatter=True,aspect='equal')
plt.grid(1)
display(plt.gcf())
#here i can remove cylinder
# wrapper around minimize.
def fit(pts,odr2,fit_func):
""" fit pts and return residuals. Info are printed. odr2 is starting guess.
typically:
fit_func=cylinder_error3 #function to be minimized
odr2=(span(xg).sum()/2.,220.,0,0.) #use nominal value for guess direction
deltaR=fit_func(pts,odr2,fit_func).
pts and deltaR enter and exit with z in microns
"""
pts=pts/[1,1,1000.] #data are um, convert to mm before fit
result=minimize(fit_func,x0=(odr2),
args=(pts,),
options={'maxiter':500},
method='Nelder-Mead')#,callback=p)
print '-----------------------------------'
print 'Results of fit on region:'
print result
odr=result.x
fom,deltaR,coeff=fit_func(odr,pts,extra=True)
deltaR[:,2]=deltaR[:,2]*1000
print 'Angle of cyl axis with y axis (deg):'
print np.arccos(np.sqrt(1-(odr[-2:]**2).sum()))*180/np.pi
print 'Radius:'
print coeff
print 'rms residuals=%s um'%(fom*1e3)
return deltaR
odr2=(span(pts[:,0]).sum()/2.,220.,0,0.) #use nominal value for guess direction
deltaR=fit(pts,odr2,fit_func=cylinder_error3) #cylinder_error3 #function to be minimized
sigma=np.std(deltaR[:,2])
print sigma
plt.clf()
plot_points(deltaR,scatter=True,aspect='equal',vmin=-sigma,vmax=sigma)
display(plt.gcf())
plt.grid(1)
#remove outliers based on deltaR from pts
sigma=np.std(deltaR[:,2])
mean=np.mean(deltaR[:,2])
mask=(deltaR[:,2]-mean)**2<(nsigma*sigma)**2 #points to keep
pts=pts[mask,:]
print sigma,np.std(deltaR[mask,2])
odr2=(span(pts[:,0]).sum()/2.,220.,0,0.) #use nominal value for guess direction
deltaR=fit(pts,odr2,fit_func=cylinder_error3) #cylinder_error3 #function to be minimized
sigma=np.std(deltaR[:,2])
print sigma
plt.clf()
plot_points(deltaR,scatter=True,aspect='equal',vmin=-sigma,vmax=sigma)
display(plt.gcf())
plt.grid(1)
f,ap1,ap2=correctedPSD(deltaR,crop=crop,step=step,nsigma=nsigma,outname=df)
plt.savefig(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.png')))
np.savetxt(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.dat')),
np.vstack([f[:],ap1[:-1],ap2[:-1]]).T)
df1=r'../OP1surface/measure_data/calibration_m3/2013_12_23/OAB150x150_01_Height.txt'
df2=r'C:\Users\Vincenzo\Google Drive\POOL\OP1surface\OAB1\xysurface150x150_Height_transformed_deltaR_4inch.dat'
df3=r'..\OP1surface\OP1\09_op1_yxsurf_Height_transformed.dat'
f1,p1,p1c=np.genfromtxt(os.path.join(outfolder,fn_add_subfix(
os.path.basename(df1),'_PSD','.dat')),
unpack=True)
f2,p2,p2c=np.genfromtxt(os.path.join(outfolder,fn_add_subfix(
os.path.basename(df2),'_PSD','.dat')),
unpack=True)
f3,p3,p3c=np.genfromtxt(os.path.join(outfolder,fn_add_subfix(
os.path.basename(df3),'_PSD','.dat')),
unpack=True)
plt.figure('PSD')
plt.clf()
plt.title('PSD of mandrels')
plt.plot(f1[1:],p1[1:],'b--',label='Ref. flat Raw')
plt.plot(f1[1:],p1c[1:],'b',label='Ref. flat outlier removed')
plt.plot(f2[1:],p2[1:],'r--',label='OAB1 Raw')
plt.plot(f2[1:],p2c[1:],'r',label='OAB1 outlier removed')
plt.plot(f3[1:],p3[1:],'g--',label='OP1 Raw')
plt.plot(f3[1:],p3c[1:],'g',label='OP1 outlier removed')
plt.loglog()
plt.ylabel('um^2')
plt.xlabel('mm^-1')
plt.xlim((0.8e-2,2e0))
plt.legend(loc=0)
plt.grid(1)
display(plt.gcf())
plt.savefig(os.path.join(outfolder,'mandrels_PSD.png'))
#corrected only
plt.figure('PSD')
plt.clf()
plt.title('PSD of mandrels')
#plt.plot(f1[1:],p1[1:],'b--',label='Ref. flat Raw')
plt.plot(f1[1:],p1c[1:],'b',label='Ref. flat')
#plt.plot(f2[1:],p2[1:],'r--',label='OAB1 Raw')
plt.plot(f2[1:],p2c[1:],'r',label='OAB1')
#plt.plot(f3[1:],p3[1:],'g--',label='OP1 Raw')
plt.plot(f3[1:],p3c[1:],'g',label='OP1')
plt.loglog()
plt.ylabel('um^2')
plt.xlabel('mm^-1')
plt.xlim((0.8e-2,2e0))
plt.legend(loc=0)
plt.grid(1)
display(plt.gcf())
plt.savefig(os.path.join(outfolder,'mandrels_PSD_corronly.png'))